スコアマッチングとその周辺

Author

March

Published

July 4, 2025

\[ \newcommand{\darrow}{\overset{d}{\to}} \newcommand{\parrow}{\overset{p}{\to}} \newcommand{\warrow}{\rightsquigarrow} \newcommand{\bm}[1]{\boldsymbol{#1}} \newcommand{\d}{\,\text{d}} \newcommand{\dd}[3][]{\frac{\text{d}^{#1}{#2}}{\text{d}^{#1}{#3}}} \newcommand{\pd}[3][]{\frac{\partial^{#1}{#2}}{\partial^{#1}{#3}}} \newcommand{\ang}[1]{\langle{#1}\rangle} \]

はじめに

この資料では,スコアマッチングを中心に,正規化定数に依らない統計的推測の方法を紹介する.以下のようなモデルを考える.

\[ q(x\mid\theta)=\frac{\tilde{q}(x\mid\theta)}{Z(\theta)} \]

ここで,\(Z(\theta)\) は正規化定数で,連続分布の場合は \[ Z(\theta)=\int\tilde{q}(x\mid\theta)\d{x} \]

として定まる.離散分布の場合は, \[ Z(\theta)=\sum_{x}\tilde{q}(x\mid\theta) \]

として定まる.この正規化定数の評価が困難な場合を考える. そのような状況では最尤推定は利用できないため,何か別の方法が必要になる. はじめに,正規化定数の評価が困難なモデルの例を紹介する.

Example 1 (エネルギーベースモデル) エネルギーベースモデルは,エネルギー関数 \(E(x\mid\theta)\) を用いて, \[ q(x\mid\theta)\propto e^{-E(x\mid\theta)} \]

と表されるモデルである.\(-\log{q}(x\mid\theta)\)\(x\) に依る項のみを \(E(x\mid\theta)\) によって直接モデリングする.機械学習分野では \(E(x\mid\theta)\) をニューラルネットなどで構成し,複雑なデータ分布を学習する.一般に正規化定数は評価できない.

Example 2 (Boltzmann マシン) Bolzmann マシンはエネルギー関数 \[ E(\bm{x}\mid \bm\theta)=\sum_{j=1}^{d}\theta_{j}x_{j}+\sum_{j<k}\theta_{jk}x_{j}x_{k} \]

による \(\mathfrak{X}=\{-1,1\}^{d}\) 上のエネルギーベースモデルである.正規化定数は \[ Z(\bm\theta)=\sum_{\bm{x}\in\{-1,1\}^{d}}e^{-E(\bm{x}\mid\bm\theta)} \]

\(2^d\) 項の和となるため,\(d\) が大きいと評価は現実的でない.

Example 3 (切断分布) \(\mathfrak{X}\) 上のモデル \(p(x\mid\theta)\) を用いて,\(V\subset\mathfrak{X}\) 上のモデルを \[ q(x\mid\theta)\propto p(x\mid\theta)1_{V}(x) \]

によって構成する.このようなモデル(分布)を切断分布という1.正規化定数 \(Z(\theta)\)\[ Z(\theta)=\int_{V}p(x\mid\theta)\d{x} \]

によって与えられる.\(p\) が正規分布の場合ですら,\(Z(\theta)\) は陽に表せない. 例えば,Liu, Kanamori, and Williams (2022) では,シカゴの犯罪発生地点のデータに多変量混合正規分布の切断分布を当てはめることを検討している.ここで,\(V\subset\mathbb{R}^{2}\) は多角形で,シカゴの形状を表現する.他にも,生存時間解析では切断分布がよく現れる.

Example 4 (切断 Gaussian グラフィカルモデル) \([0,\infty)^{d}\) 上の切断 Gaussian グラフィカルモデルは \[ q(\bm{x}\mid\Sigma)\propto \exp\left(-\frac{1}{2}\bm{x}\Sigma^{-1}\bm{x}\right) \]

で表されるモデルである.\(\Sigma\) を推定することで,変数間の依存関係を調べることができる. Lin, Drton, and Shojaie (2016) では,RNA-seq データにこのモデルを適用した.RNA-seq データは数百種類の遺伝子の発現量を並べた非負ベクトルデータである.複雑な疾病を解析するには,遺伝子発現の相互作用を調べることが重要(らしい).

Example 5 (トーラスグラフモデル) \(d\) 次元トーラスは円周 \(\mathbb{S}^1\)\(d\) 個の直積である.角度の \(d\) 変量データ \(\bm{x}\in[0,2\pi)^{d}\)\(d\) 次元トーラス上の1点とみなせる.

Klein et al. (2020) によって導入されたトーラスグラフモデルは \[ q(\bm{x}\mid\bm{\theta})\propto \exp\left( \sum_{j=1}^{d}\begin{bmatrix}\theta_{j,1}\\\theta_{j,2}\end{bmatrix}^{\top}\begin{bmatrix}\cos{x_j}\\\sin{x_j}\end{bmatrix} +\sum_{j<k} \begin{bmatrix}\theta_{jk,1}\\\theta_{jk,2}\\\theta_{jk,3}\\\theta_{jk,4}\end{bmatrix}^\top \begin{bmatrix}\cos(x_j-x_k)\\\sin(x_j-x_k)\\\cos(x_j+x_k)\\\sin(x_j+x_k)\end{bmatrix} \right) \]

で表されるモデルで,トーラス上のグラフィカルモデルである.正規化定数の評価は見るからに困難である.

Example 6 (多様体上の確率分布) 多様体上の確率分布の多くは,評価困難な正規化定数をもつ. 例えば,超球面 \(\mathbb{S}^{d-1}\) 上の Fisher-Bingham 分布は密度 \[ q(\bm{x}\mid\bm\mu,A)\propto \exp\left( \bm{x}^\top \bm{\mu} + \bm{x}^\top A\bm{x} \right) \]

をもつ分布である.\(\mathbb{R}^d\) 上の多変量正規分布を球面上に制限することで得られる.

スコアマッチング

スコアマッチングは Hyvärinen (2005) によって提案された正規化定数に依らない統計的推測の方法である.

Fisher ダイバージェンス

Definition 1 (Fisher ダイバージェンス) \(\mathfrak{X}\subset\mathbb{R}^{d}\) 上の2つの確率分布 \(P,Q\) は,微分可能な正の密度 \(p(\bm{x}),q(\bm{x})\) をもつとする. \[ J(P\,\|\,Q):=\int_{\mathfrak{X}} p(\bm{x})\|\nabla_{\bm{x}}\log{p}(\bm{x})-\nabla_{\bm{x}}\log{q}(\bm{x})\|^{2}\d{\bm{x}} \]

を Fisher ダイバージェンスという.

Proposition 1 \(\mathfrak{X}\) が開連結なら,Fisher ダイバージェンスはダイバージェンスの性質を満たす.すなわち,

  • \(J(P\,\|\,Q)\geq0\)
  • \(J(P\,\|\,Q)=0\iff P=Q\)

が成立.

証明

Proof. 非負性と \(P=Q\Rightarrow J(P\,\|\,Q)=0\) は明らか.\(J(P\,\|\,Q)=0\) を仮定する. このとき a.e で, \(\|\nabla_{\bm{x}}\log{p}(\bm{x})-\nabla_{\bm{x}}\log{q}(\bm{x})\|^{2}=0\),すなわち \(\nabla_{\bm{x}}\log{p}(\boldsymbol{x})=\nabla_{\bm{x}}\log{q}(\bm{x})\) が成立.このとき,ある定数 \(C\) が存在して a.e で \[ \log{p}(\bm{x})=\log{q}(\bm{x})+C \]

が成り立つが,\(p,q\) は密度関数なので,\(C=0\) となる.よって,a.e で \(p=q\) が成立し,\(P=Q\)

Example 7 (\(J(P\,\|\,Q)=0\Rightarrow P=Q\) が成立しない例) 各連結成分で, \(\log p(\bm{x})-\log q(\bm{x})\) が異なる定数になるようにすることで反例が構成できる.\(\mathfrak{X}\) は開だが連結でない,すなわち,2つの非交差な開集合 \(G_{1},G_{2}\) によって, \(\mathfrak{X}=G_{1}\sqcup G_{2}\) と表せると仮定する.ここで,微分可能な \(f_1,f_2\) を用いて,\(P,Q\) の密度をそれぞれ \[ \begin{align*} p(\bm{x})&:=\alpha f_{1}(\bm{x})1_{G_{1}}(\bm{x})+(1-\alpha)f_{2}(\bm{x})1_{G_{2}}(\bm{x})\\ q(\bm{x})&:=\beta f_{1}(\bm{x})1_{G_{1}}(\bm{x})+(1-\beta)f_{2}(\bm{x})1_{G_{2}}(\bm{x}) \end{align*} \]

で定めると,\(J(P\,\|\,Q)=0\) が任意の \(\alpha,\beta\in[0,1]\) で成立.

\(\nabla_{\bm{x}}\log{q}(\bm{x})\) を(データ \(\bm{x}\) に関する)スコア関数という.スコアマッチングとは,Fisher ダイバージェンスの最小化によって(スコア関数間の距離を近づけることによって),\(P\)\(Q\) で近似する方法である.

いま,\(P\) は未知として,標本 \(\bm{X}_{1},\dots,\bm{X}_{n}\sim P\) による経験分布 \(\mathbb{P}_{n}=n^{-1}\sum_{i=1}^{n}\delta_{\bm{X}_{i}}\) との Fisher ダイバージェンスを最小化するような \(Q\) によって \(P\) を近似することを考えると,\(\nabla_{\bm{x}}\log{p}(\bm{x})\) を知ることができないため困ってしまう.次節で導出する Hyvärinen スコアは適当な条件の下で,\(\nabla_{\bm{x}}\log{p}\) の評価を回避して経験 Fisher ダイバージェンスを最小化する方法を与える.

Hyvärinen スコアの導出

Hyvärinen スコアの導出2には Green の定理(あるいは部分積分の公式)

\[ \int_{\mathfrak{X}}\partial_{j}(f(\bm{x})g(\bm{x}))\d{\bm{x}}=\int_{\partial\mathfrak{X}}f(\bm{x})g(\bm{x})\nu_{j}(\bm{x})\d{s} \]

が本質的な役割を果たす3.ただし,\(\partial_j=\partial/\partial{x_j}\) で,\((\nu_1,\dots,\nu_d)^\top\)\(\partial{\mathfrak{X}}\) の単位法線ベクトル.

Fisher ダイバージェンスは以下のように展開できる. \[ \begin{align*} J(P\,\|\,Q) &=\int_{\mathfrak{X}} p(\bm{x})\|\nabla_{\bm{x}}\log{p}(\bm{x})-\nabla_{\bm{x}}\log{q}(\bm{x})\|^{2}\d{\bm{x}}\\ &=\int_{\mathfrak{X}} p(\bm{x}) \left( \|\nabla_{\bm{x}}\log{p}(\bm{x})\|^{2} +\|\nabla_{\bm{x}}\log{q}(\bm{x})\|^{2} -2\langle\nabla_{\bm{x}}\log{p}(\bm{x}),\nabla_{\bm{x}}\log{q}(\bm{x})\rangle \right)\d{\bm{x}} \end{align*} \]

ここで,第3項は Green の定理によって

\[ \begin{align*} &\int_\mathfrak{X}p(\bm{x})\langle\nabla_{\bm{x}}\log{p}(\bm{x}),\nabla_{\bm{x}}\log{q}(\bm{x})\rangle\d{x}\\ &=\sum\limits_{j=1}^{d}\int_\mathfrak{X}p(\bm{x})\frac{\partial\log{p}(\bm{x})}{\partial{x_{j}}}\frac{\partial\log{q}(\bm{x})}{\partial{x_{j}}}\d{x}\\ &=\sum\limits_{j=1}^{d}\int_\mathfrak{X}\frac{\partial{p}(\bm{x})}{\partial{x_{j}}}\frac{\partial\log{q}(\bm{x})}{\partial{x_{j}}}\d{x}\\ &=-\sum\limits_{j=1}^{d}\int_\mathfrak{X}p(\bm{x})\frac{\partial^{2}\log{q}(\bm{x})}{\partial{x_{j}^{2}}}\d{x}+\sum\limits_{j=1}^{d}\int_{\partial{\mathfrak{X}}}p(\bm{x})\frac{\partial\log{q}(\bm{x})}{\partial{x_{j}}}\nu_{j}(\bm{x})\d{s} \end{align*} \]

と展開できる4.このとき,境界項が無視できる(0になる)ならば,

\[ \begin{align*} J(P\,\|\,Q) &=E_{P}\Big[\|\nabla_{\bm{x}}\log{q\|^{2}+2\Delta_{\bm{x}}\log{q}}\Big]+E_{P}\Big[\|\nabla_{\bm{x}}\log{p}\|^{2}\Big] \end{align*} \]

となる.ただし,\(E_P[\,\cdot\,]\)\(P\) による期待値を意味する.\(\Delta_{\bm{x}}\) は Laplacian で,\(\Delta_{\bm{x}}=\sum_{j=1}^{d}\partial^{2}/\partial{x_{j}^{2}}\) である.\(\|\nabla_{\bm{x}}\log{q\|^{2}+2\Delta_{\bm{x}}\log{q}}\) を Hyvärinen スコアという.

経験分布 \(\mathbb{P}_{n}\) との Fisher ダイバージェンスを最小化するような \(Q\) によって \(P\) を近似することを考える.\(E_{P}[\|\nabla_{X}\log{p}\|^{2}]\)\(Q\) と無関係なので,Fisher ダイバージェンス の最小化は Hyvärinen スコアの期待値の最小化である.したがって,経験 Fisher ダイバージェンスの最小化は \[ E_{\mathbb{P}_{n}}\Big[\|\nabla_{\bm{x}}\log{q\|^{2}+2\Delta_{\bm{x}}\log{q}}\Big]= \sum_{i=1}^{n}\left(\|\nabla_{x}\log{q}(\bm{X}_{i})\|^{2}+2\Delta_{\bm{x}}\log{q}(\bm{X}_i)\right) \]

の最小化によって達成される.

経験 Fisher ダイバージェンスを最小化する推定量をスコアマッチング推定量という.スコアマッチング推定量の漸近的な振る舞いは M 推定量の枠組みで議論される.

Example 8 \(\mathfrak{X}=(a_{1},b_{1})\times\cdots\times(a_{d},b_{d})\) の場合,境界項を無視できる条件は, \[ \lim_{x_j\searrow a_j}p(\bm{x})\frac{\partial\log{q}(\bm{x})}{\partial{x_j}}=\lim_{x_j\nearrow b_j}p(\bm{x})\frac{\partial\log{q}(\bm{x})}{\partial{x_j}} \]

が(両辺有限で)任意の \(j=1,\dots,d\) で成立することである.

Example 9 (境界項が無視できる例 トーラスグラフモデル) \(p\) が連続なら,トーラスグラフモデルの境界項は0になる.

Example 10 (境界項が無視できない例 その1 非負分布) \(P=\text{Exp}(\theta_0),Q=\text{Exp}(\theta)\) とする.このとき, \[ \begin{align*} \lim_{x\searrow0}p(x)\dd{\log{q}(x)}{x} &=\lim_{x\searrow0}\frac{1}{\theta_0}e^{-x/\theta_0}\dd{}{x}(-\frac{x}{\theta}-\log\theta) =\lim_{x\searrow0}-\frac{1}{\theta\theta_0}e^{-x/\theta_0} =-\frac{1}{\theta\theta_0}\\ \lim_{x\nearrow\infty}p(x)\dd{\log{q}(x)}{x} &=\lim_{x\nearrow\infty}-\frac{1}{\theta\theta_0}e^{-x/\theta_0} =0 \end{align*} \]

となり,境界項は消えない.

Example 11 (境界項が無視できない例 その2 切断分布) \(P\)\((a,b)\) 上に制限された正規分布 \(N(\mu_0,\sigma_0^2)\)\(Q\) も同様に制限された正規分布 \(N(\mu,\sigma^2)\) とする.密度は \[ q(x\mid\mu,\sigma^2)=\frac{1}{C(\mu,\sigma^2)}\frac{1}{\sqrt{2\pi\sigma^2}}\exp\left\{-\frac{(x-\mu)^2}{2\sigma^2}\right\},\quad x\in(a,b) \]

である.ただし, \[ C(\mu,\sigma^2)=\int_{a}^{b}\frac{1}{\sqrt{2\pi\sigma^2}}\exp\left\{-\frac{(x-\mu)^2}{2\sigma^2}\right\}\d{x} \]

である.このとき, \[ \dd{\log{q}(x)}{x}=\dd{}{x}-\frac{(x-\mu)^2}{2\sigma^2}=-\frac{x-\mu}{\sigma^2} \]

となり,境界項が消えるとは限らない.

Example 12 (境界項が無視できない例 その3 裾が重い分布と裾が軽い分布) \(P=\text{Cauchy}(0,1)\) で,\(Q\) は密度 \(q(x)\propto\exp(-(x-\mu)^4/4)\) をもつとする. \[ \dd{\log{q}(x)}{x}=-(x-\mu)^3 \]

より, \[ \lim_{x\nearrow\infty}p(x)\dd{\log{q}(x)}{x}= \lim_{x\nearrow\infty}-\frac{1}{\pi}\frac{(x-\mu)^3}{1+x^2}=\infty \]

となり,境界項が発散してしまう.

Example 13 (スコアマッチング推定量の例 その1 正規分布) 多変量正規分布 \(N(\bm{\mu},\Sigma)\) のスコアマッチング推定量を求める.

途中式

\(\log q(\bm{x}\mid\bm{\mu},\Sigma)\) の勾配と Laplacian はそれぞれ \[ \begin{align*} \nabla_{\bm{x}}\log{q}(\bm{x}\mid\bm{\mu},\Sigma) &=\nabla_{\bm{x}}\left\{-\frac{1}{2}(\bm{x}-\bm{\mu})^{\top}\Sigma^{-1}(\bm{x}-\bm{\mu})\right\} =-\Sigma^{-1}(\bm{x}-\bm{\mu})\\[5pt] \Delta_{\bm{x}}\log{q}(\bm{x}\mid\bm{\mu},\Sigma)&=-\mathrm{tr}\,\Sigma^{-1} \end{align*} \]

と表される.したがって,Hyvärinen スコアは \[ H(\bm{x}\mid\bm{\mu},\Sigma)-2\mathrm{tr}\,\Sigma^{-1}+(\bm{x}-\bm{\mu})^{\top}\Sigma^{-2}(\bm{x}-\bm{\mu}) \]

である.\(H(\bm{x}\mid\bm{\mu},\Sigma)\) の勾配を求める. \[ \begin{align*} \nabla_{\bm{\mu}}H(\bm{x}\mid\bm{\mu},\Sigma) &=\nabla_{\bm{\mu}}\left((\bm{x}-\bm{\mu})^{\top}\Sigma^{-2}(\bm{x}-\bm{\mu})\right)\\ &=2\Sigma^{-2}(\bm{\mu}-\bm{x})\\[5pt] \nabla_{\Sigma^{-1}}H(\bm{x}\mid\bm{\mu},\Sigma) &=-2\nabla_{\Sigma^{-1}}\mathrm{tr}\,\Sigma^{-1}+\nabla_{\Sigma^{-1}}(\bm{x}-\bm{\mu})^{\top}\Sigma^{-2}(\bm{x}-\bm{\mu})\\ &=-2I_{d}+2(\bm{x}-\bm{\mu})(\bm{x}-\bm{\mu})^{\top}\Sigma^{-1} \end{align*} \]

ただし,\(\nabla_{A}\bm{x}^{\top}A^{2}\bm{x}=\bm{x}\bm{x}^{\top}A^{\top}+A^{\top}\bm{x}\bm{x}^{\top}\) を用いた.

標本 \(\bm{X}_1,\dots,\bm{X}_n\) に対し,推定方程式 \[ \begin{align*} \sum_{i=1}^{n}\nabla_{\bm{\mu}}H(\bm{X}_i\mid\bm{\mu},\Sigma)&=\bm{0}\\ \sum_{i=1}^{n}\nabla_{\Sigma^{-1}}H(\bm{X}_i\mid\bm{\mu},\Sigma)&=O_d \end{align*} \]

を解いてスコアマッチング推定量を求める.

\[ \hat{\bm\mu}_{n}=\frac{1}{n}\sum_{i=1}^{n}\bm{X}_{i},\enspace \hat\Sigma_n=\frac{1}{n}\sum_{i=1}^{n}(\bm{X}_{i}-\hat{\bm{\mu}}_{n})(\bm{X}_{i}-\hat{\bm{\mu}}_{n})^{\top} \]

となり,最尤推定量と一致する.

Example 14 (スコアマッチング推定量の例 その2 指数型分布族) 以下の密度 \(q(\bm{x}\mid\bm{\theta})\) をもつモデルのスコアマッチング推定量を求める.ただし,\(b:\mathbb{R}^{d}\to\mathbb{R},\bm{t}:\mathbb{R}^{d}\to\mathbb{R}^{m}\) である. \[ q(\bm{x}\mid\bm{\theta})=\frac{1}{Z(\bm{\theta})}b(\bm{x})e^{-\bm{\theta}^{\top}\bm{t}(\bm{x})} \]

途中式

\[ \begin{align*} B_{1}(\bm{x})&=(\partial_{1}\log{b}(\bm{x}),\dots,\partial_{d}\log{b}(\bm{x}))^{\top}\\ B_{2}(\bm{x})&=(\partial_{1}^{2}\log{b}(\bm{x}),\dots,\partial_{d}^{2}\log{b}(\bm{x}))^{\top}\\ T_{1}(\bm{x})&=\begin{bmatrix} \pd{t_{1}}{x_{1}}&\cdots&\pd{t_{1}}{x_{d}}\\ \vdots&\ddots&\vdots\\ \pd{t_{m}}{x_{1}}&\cdots&\pd{t_{m}}{x_{d}} \end{bmatrix}\\ T_{2}(\bm{x})&=\begin{bmatrix} \pd[2]{t_{1}}{x_{1}}&\cdots&\pd[2]{t_{1}}{x_{d}}\\ \vdots&\ddots&\vdots\\ \pd[2]{t_{m}}{x_{1}}&\cdots&\pd[2]{t_{m}}{x_{d}} \end{bmatrix} \end{align*} \]

とすると, \[ \begin{align*} \nabla_{\bm{x}}\log{q}(\bm{x}) &=B_{1}(\bm{x})-\bm{\theta}^{\top}T_{1}(\bm{x})\\ \Delta_{\bm{x}}\log{q}(\bm{x}) &=B_{2}(\bm{x})^\top\bm{1}_{d}-\bm{\theta}^{\top}T_{2}(\bm{x})\bm{1}_{d} \end{align*} \]

より,Hyvärinen スコアは \[ H(\bm{x}\mid\bm{\theta})=B_{1}(\bm{x})^{\top}B_{1}(\bm{x})+\bm{\theta}^{\top}T_{1}(\bm{x})T_{1}(\bm{x})^{\top}\bm{\theta}-2\bm{\theta}^{\top}T_{1}(\bm{x})B_{1}(\bm{x})+2B_{2}(\bm{x})^{\top}\bm{1}_{d}-2\bm{\theta}^{\top}T_{2}(\bm{x})\bm{1}_{d} \]

となる.\(\nabla_{\bm{\theta}}H(\bm{x}\mid\bm{\theta})\)\[ \nabla_{\bm{\theta}}H(\bm{x}\mid\bm{\theta}) =2T_{1}(\bm{x})T_{1}(\bm{x})^{\top}\bm{\theta}-2T_{1}(\bm{x})B_{1}(\bm{x})-2T_{2}(\bm{x})\bm{1}_{d} \]

となる.したがって,

標本 \(\bm{X}_{1},\dots,\bm{X}_{n}\) に関するスコアマッチング推定量は \[ \hat{\bm{\theta}}_{n} =\left(\sum_{i=1}^{n}T_{1}(\bm{X}_{i})T_{1}(\bm{X}_{i})^{\top}\right)^{-1} \left\{\sum_{i=1}^{n}\left(T_{1}(\bm{X}_{i})B_{1}(\bm{X}_{i})+T_{2}(\bm{X}_{i})\bm{1}_{d}\right)\right\} \]

と,線型方程式の解として定まる.

生成モデルへの応用

スコアマッチングの主な応用先の1つとして,生成モデルがある.

Langevin モンテカルロ

密度 \(p(\bm{x})\) をもつ分布 \(P\) からデータを生成したい.Langevin モンテカルロは,適当な初期値 \(\bm{x}_{0}\) から始めて更新式 \[ \bm{x}_{t+1}:=\bm{x}_{t}+\varepsilon\nabla_{\bm{x}}\log{p}(\bm{x}_{t})+\sqrt{2\varepsilon}\bm{z}_{t},\quad \bm{z}_{t}\sim N(0,I) \]

によって \(\bm{x}_{T}\) を得るアルゴリズムである.\(T\to\infty,\varepsilon\to0\) の極限で,\(\bm{x}_{T}\)\(P\) に従うことが知られている.

生成モデルとスコアマッチング

Langevin モンテカルロの更新式から,\(P\) からデータを生成するためには,\(\nabla_{\bm{x}}\log{p}(\bm{x})\) を学習すればよいことがわかる. 損失関数として, \[ E_{P}\Big[\|\bm\psi(\bm{x})-\nabla_{\bm{x}}\log{p}(\bm{x})\|^{2}\Big] \]

を用いると,これは Fisher ダイバージェンスであり,正則条件の下で,Hyvärinen スコア \[ E_{P}\Big[\|\bm\psi(\bm{x})\|^{2}+2\operatorname{tr}\nabla_{\bm{x}}\bm\psi(\bm{x})^{\top}\Big] \]

の最小化と等しくなる.

Hyvärinen スコアの課題

画像データのように,\(\bm{x}\) の次元 \(d\) が大きいと,\(\operatorname{tr}\nabla_{\bm{x}}\bm\psi^{\top}\) の数値計算が大変になる.

関数 \(\bm{f}:\mathbb{R}^{l}\to\mathbb{R}^{m}\) の自動微分を考えよう.\((y_{1},\dots,y_{m})=\bm{f}(x_{1},\dots,x_{l})\) とする.前向きの自動微分では,Jacobi 行列 \[ \begin{bmatrix} \frac{\partial y_{1}}{\partial x_{1}}&\cdots&\frac{\partial{y_{1}}}{\partial{x_{l}}}\\ \vdots&\ddots&\vdots\\ \frac{\partial{y}_{m}}{\partial{x}_{1}}&\dots&\frac{\partial{y}_{m}}{\partial{x}_{l}} \end{bmatrix} \]

の任意の1列を一度に計算できる.後向きの自動微分では,Jacobi 行列の任意の1行を一度に計算できる.\(\operatorname{tr}\nabla_{\bm{x}}\bm\psi^{\top}\) を求めるには,Jacobi 行列の対角成分が必要で,\(d\) 回の自動微分が要求される.とても効率が悪い.

デノイジング・スコアマッチング

Vincent (2011) によって提案されたデノイジング・スコアマッチングは \(\operatorname{tr}\bm\nabla_{\bm{x}}\bm\psi^{\top}\) の計算を回避できる方法の1つである.いま,\(\bm{x}\) は興味ある分布 \(P\) に従うとして,\(\tilde{\bm{x}}=\bm{x}+\bm{z}\) とする.\(\bm{z}\) は既知の分布に従うノイズである.

Warning

この章では,\(p(\bm{u})\) は引数 \(\bm{u}\) が従う分布の密度関数として,\(E_{p(\bm{u})}[\,\cdot\,]\)\(\bm{u}\) が従う分布に関する期待値を意味する.

\(\tilde{\bm{x}}\) に関するスコアマッチングを考える(\(\bm{z}\) が十分小さいノイズなら \(\bm{x}\) に関するスコアマッチングと思える).スコアのモデル \(\bm\psi\) に関する Fisher ダイバージェンスは \[ E_{p(\tilde{\bm{x}})}\Big[\|\bm\psi(\tilde{\bm{x}})-\nabla_{\tilde{\bm{x}}}\log{p}(\tilde{\bm{x}})\|^{2}\Big] \]

で与えられる.Hyvärinen スコアは \[ E_{p(\tilde{\bm{x}})}\big[\|\bm\psi(\tilde{\bm{x}})\|^{2}+2\operatorname{tr}\nabla_{\tilde{\bm{x}}}\bm\psi(\tilde{\bm{x}})^{\top}\big] \]

だが,ここでは別の表現を与える.

途中式

Fisher ダイバージェンスを展開すると, \[ \begin{align*} &E_{p(\tilde{\bm{x}})}\Big[\|\bm\psi(\tilde{\bm{x}})-\nabla_{\tilde{\bm{x}}}\log{p}(\tilde{\bm{x}})\|^{2}\Big]\\ &=E_{p(\tilde{\bm{x}})}\Big[\|\bm\psi(\tilde{\bm{x}})\|^{2}\Big]+E_{p(\tilde{\bm{x}})}\Big[\|\nabla_{\tilde{\bm{x}}}\log{p}(\tilde{\bm{x}})\|^{2}\Big]-2E_{p(\tilde{\bm{x}})}\Big[\big\langle\bm\psi(\tilde{\bm{x}}),\nabla_{\tilde{\bm{x}}}\log{p}(\tilde{\bm{x}})\big\rangle\Big] \end{align*} \]

となる.第3項を評価する. \[ \begin{align*} \nabla_{\tilde{\bm{x}}}\,p(\tilde{\bm{x}}) &=\nabla_{\tilde{\bm{x}}}\int p(\tilde{\bm{x}}\mid\bm{x})p(\bm{x})\d{\bm{x}}\\ &=\int \nabla_{\tilde{\bm{x}}}\,p(\tilde{\bm{x}}\mid\bm{x})p(\bm{x})\d{\bm{x}}\\ &=\int \nabla_{\tilde{\bm{x}}}\log{p}(\tilde{\bm{x}}\mid\bm{x})p(\tilde{\bm{x}}\mid\bm{x})p(\bm{x})\d{\bm{x}} \end{align*} \]

を認めると, \[ \begin{align*} \int\big\langle\bm\psi(\tilde{\bm{x}}),\nabla_{\tilde{\bm{x}}}\log{p}(\tilde{\bm{x}})\big\rangle p(\tilde{\bm{x}})\d\tilde{\bm{x}} &=\int\big\langle\bm\psi(\tilde{\bm{x}}),\nabla_{\tilde{\bm{x}}}p(\tilde{\bm{x}})\big\rangle\d\tilde{\bm{x}}\\ &=\int\int\big\langle\bm\psi(\tilde{\bm{x}}),\nabla_{\tilde{\bm{x}}}\log{p}(\tilde{\bm{x}}\mid \bm{x})\big\rangle p(\tilde{\bm{x}}\mid\bm{x})p(\bm{x})\d{\bm{x}}\d\tilde{\bm{x}}\\ &=E_{p(\tilde{\bm{x}}\mid\bm{x})p(\bm{x})}\Big[\big\langle\bm\psi(\tilde{\bm{x}}),\nabla_{\tilde{\bm{x}}}\log{p}(\tilde{\bm{x}}\mid\bm{x})\big\rangle\Big] \end{align*} \]

が得られる.したがって, \[ \begin{align*} &E_{p(\tilde{\bm{x}})}\Big[\|\bm\psi(\tilde{\bm{x}})-\nabla_{\tilde{\bm{x}}}\log{p}(\tilde{\bm{x}})\|^{2}\Big]\\ &=E_{p(\tilde{\bm{x}})}\Big[\|\bm\psi(\tilde{\bm{x}})\|^{2}\Big]+E_{p(\tilde{\bm{x}})}\Big[\|\bm\psi(\tilde{\bm{x}})\|^{2}\Big]-2E_{p(\tilde{\bm{x}})}\Big[\big\langle\bm\psi(\tilde{\bm{x}}),\nabla_{\tilde{\bm{x}}}\log{p}(\tilde{\bm{x}})\big\rangle\Big]\\ &=E_{p(\tilde{\bm{x}}|\bm{x})p(\bm{x})}\Big[\|\bm\psi(\tilde{\bm{x}})\|^{2}\Big]+E_{p(\tilde{\bm{x}})}\Big[\|\nabla_{\tilde{\bm{x}}}\log{p}(\tilde{\bm{x}})\|^{2}\Big]-2E_{p(\tilde{\bm{x}}\mid\bm{x})p(\bm{x})}\Big[\big\langle\bm\psi(\tilde{\bm{x}}),\nabla_{\tilde{\bm{x}}}\log{p}(\tilde{\bm{x}}\mid\bm{x})\big\rangle\Big]\\ &=E_{p(\tilde{\bm{x}}|\bm{x})p(\bm{x})}\Big[\|\bm\psi(\tilde{\bm{x}})-\nabla_{\tilde{\bm{x}}}\log{p}(\tilde{\bm{x}}\mid\bm{x})\|^{2}\Big]\\ &\hphantom{=}+E_{p(\tilde{\bm{x}})}\Big[\|\nabla_{\tilde{\bm{x}}}\log{p}(\tilde{\bm{x}})\|^{2}\Big]-E_{p(\tilde{\bm{x}}\mid\bm{x})p(\bm{x})}\Big[\|\nabla_{\tilde{\bm{x}}}\log{p}(\tilde{\bm{x}}\mid\bm{x})\|^{2}\Big] \end{align*} \]

となる.第1項以外は,\(\bm\psi\) に依らないので無視できる.よって,

\(E_{p(\tilde{\bm{x}})}\Big[\|\bm\psi(\tilde{\bm{x}})-\nabla_{\tilde{\bm{x}}}\log{p}(\tilde{\bm{x}})\|^{2}\Big]\) の最小化は \(E_{p(\tilde{\bm{x}}|\bm{x})p(\bm{x})}\Big[\|\bm\psi(\tilde{\bm{x}})-\nabla_{\tilde{\bm{x}}}\log{p}(\tilde{\bm{x}}\mid\bm{x})\|^{2}\Big]\) の最小化と等価になる. \(\nabla_{\tilde{\bm{x}}}\log p(\tilde{\bm{x}}\mid\bm{x})\) が扱いやすいように,ノイズ \(\bm{z}\) を設計すると,後者の最小化は Hyvärinen スコアの最小化に比べて容易になる.

デノイジング・オートエンコーダとデノイジング・スコアマッチング

デノイジング・オートエンコーダ(DAE)はエンコーダ \(\bm{f}_{e}\) とデコーダ \(\bm{f}_{d}\) を期待損失 \[ E_{p(\tilde{\bm{x}}|\bm{x})p(\bm{x})}\Big[\|\bm{f}_{d}\circ\bm{f}_{e}(\tilde{\bm{x}})-\bm{x}\|^{2}\Big] \]

を最小化するように学習する. つまり,ノイズをよく除去するような \(\bm{f}_{e},\bm{f}_{d}\) を学習する.いま,デノイジング・スコアマッチングにおいて \(\bm{z}\sim N(0,\sigma^{2}I)\) とすると,\(\nabla_{\tilde{\bm{x}}}\log{p}(\tilde{\bm{x}}\mid\bm{x})=\nabla_{\tilde{\bm{x}}}(-\|\tilde{\bm{x}}-\bm{x}\|^{2}/2\sigma^{2})=(\bm{x}-\tilde{\bm{x}})/\sigma^{2}\) であり,期待損失は \[ E_{p(\tilde{\bm{x}}|\bm{x})p(\bm{x})}\Big[\|\bm\psi(\tilde{\bm{x}})-\nabla_{\tilde{\bm{x}}}\log{p}(\tilde{\bm{x}}\mid\bm{x})\|^{2}\Big] =\sigma^{-4}E_{p(\tilde{\bm{x}}|\bm{x})p(\bm{x})}\Big[\|\sigma^{2}\bm\psi(\tilde{\bm{x}})+\tilde{\bm{x}}-\bm{x}\|^{2}\Big] \]

となる. \[ \bm{f}_{d}\circ\bm{f}_{e}(\tilde{\bm{x}})=\sigma^{2}\bm\psi(\tilde{\bm{x}})+\tilde{\bm{x}} \]

と思うと,DAE は Fisher ダイバージェンスを最小化するようにスコアを学習していると思える.

スライス・スコアマッチング

Song et al. (2019) は,\(\operatorname{tr}\nabla_{\bm{x}}\bm\psi^{\top}\) の評価を回避するための方法として,スライス・スコアマッチングを提案した.スライス・スコアマッチングはスライス Wasserstein にインスパイアされた方法で,期待損失として \[ E_{p(\bm{v})p(\bm{x})}\Big[\big(\bm{v}^{\top}\bm\psi(\bm{x})-\bm{v}^{\top}\nabla_{\bm{x}}\log{p}(\bm{x})\big)^{2}\Big] \]

を考える.この損失も Hyvärinen スコアのように,\(p\) に依存しない形に書き換えることができる.

途中式

\[ \begin{align*} &E_{p(\bm{x})}\Big[\big(\bm{v}^{\top}\bm\psi(\bm{x})-\bm{v}^{\top}\nabla_{\bm{x}}\log{p}(\bm{x})\big)^{2}\Big]\\ &=E_{p(\bm{x})}\Big[(\bm{v}^{\top}\bm\psi(\bm{x}))^{2}\Big]+E_{p(\bm{x})}\Big[(\bm{v}^{\top}\nabla_{\bm{x}}\log{p}(\bm{x}))^{2}\Big]-2E_{p(\bm{x})}\Big[(\bm{v}^{\top}\nabla_{\bm{x}}\log{q}(\bm{x}))(\bm{v}^{\top}\nabla_{\bm{x}}\log{p}(\bm{x}))\Big]\\ \end{align*} \]

と分解し,例によって第3項を評価する. \[ \begin{align*} E_{p(\bm{x})}\Big[(\bm{v}^{\top}\bm\psi(\bm{x}))(\bm{v}^{\top}\nabla_{\bm{x}}\log{p}(\bm{x}))\Big] &=\int(\bm{v}^{\top}\bm\psi(\bm{x}))(\bm{v}^{\top}\nabla_{\bm{x}}\log{p}(\bm{x}))p(\bm{x})\d{\bm{x}}\\ &=\int(\bm{v}^{\top}\bm\psi(\bm{x}))(\bm{v}^{\top}\nabla_{\bm{x}}\,{p}(\bm{x}))\d{\bm{x}}\\ &=\sum\limits_{j=1}^{d}\int(\bm{v}^{\top}\bm\psi(\bm{x}))v_{j}\partial_{j}\,{p}(\bm{x})\d{\bm{x}}\\ &=-\sum\limits_{j=1}^{d}\int\partial_{j}(\bm{v}^{\top}\bm\psi(x))v_{j}{p}(\bm{x})\d{\bm{x}}+\text{boundary terms}\\ &=-\int(\bm{v}^{\top}\nabla_{\bm{x}}\bm\psi(x)^{\top}\bm{v})p(\bm{x})\d{\bm{x}}+\text{boundary terms}\\ &=-E_{p(\bm{x})}\Big[\bm{v}^{\top}\nabla_{\bm{x}}\bm\psi(\bm{x})^{\top}\bm{v}\Big]+\text{boundary terms} \end{align*} \]

境界項(の \(\bm{v}\) に関する期待値)が消えることを仮定すると, \[ \begin{align*} &E_{p(\bm{v})p(\bm{x})}\Big[\big(\bm{v}^{\top}\bm\psi(\bm{x})-\bm{v}^{\top}\nabla_{\bm{x}}\log{p}(\bm{x})\big)^{2}\Big]\\ &=E_{p(\bm{v})p(\bm{x})}\Big[(\bm{v}^{\top}\bm\psi(\bm{x}))^{2}\Big]+2E_{p(\bm{v})p(\bm{x})}\Big[\bm{v}^{\top}\nabla_{\bm{x}}\bm\psi(\bm{x})^{\top}\bm{v}\Big]+C \end{align*} \]

となる.

実際には経験損失 \[ \sum\limits_{m=1}^{M}\sum\limits_{i=1}^{n}\Big[(\bm{v}_{i,m}^{\top}\bm\psi(\bm{x}_{i}))^{2}+2\bm{v}_{i,m}^{\top}\nabla_{\bm{x}}\bm\psi(\bm{x}_{i})^{\top}\bm{v}_{i,m}\Big] \]

を最小化することになるが,Song et al. (2019) では \(M=1\) で十分良いと大胆な主張をしている.

Caution

多分 \(n\) が十分大きくないと \(M=1\) では機能しない.

分散補正付きスライス・スコアマッチング

Caution

本節では,通常のスコアマッチングへの適用も想定して,\(\bm\psi\)\(\nabla_{\bm{x}}\log{q}\) とおく.

Song et al. (2019) は,分散を補正した損失関数として, \[ \|\nabla_{\bm{x}}\log{q}(\bm{x})\|^{2}+2\bm{v}^{\top}\operatorname{Hess}_{\bm{x}}\log{q}(\bm{x})\bm{v} \]

を用いることも提案している(トレースの部分だけスライスする).これは Hutchinson’s trick を使った Laplacian の近似に他ならない.

Proposition 2 (Hutchinson’s trick) \(E[\bm{v}\bm{v}^{\top}]=I\) とすると,\(\operatorname{tr}A=E[\bm{v}^{\top}A\bm{v}]\) である.

証明

Proof. \[ E[\bm{v}^{\top}A\bm{v}]=\sum_{j=1}^{d}\sum_{k=1}^{d}A_{jk}E[v_{j}v_{k}]=\sum_{j=1}^{d}A_{jj}=\operatorname{tr}A. \]

Hutchinson’s trick を使って計算コストを削減できることがあるので,有用.今回のような状況でも \[ \bm{v}^{\top}\operatorname{Hess}_{\bm{x}}f(\bm{x})\bm{v}=\bm{v}^{\top}\nabla_{\bm{x}}(\bm{v}^{\top}\nabla_{\bm{x}}f(\bm{x})) \]

と表すと,自動微分の回数が2回で済むので効率が良い.

Caution

一方で,\(\bm{v}\) のモンテカルロサンプルのサイズによってはかえって非効率になるような気もする.

Langevin モンテカルロのバイアス補正

デノイジング・スコアマッチングで学習したスコアにはバイアスがある.デノイジング・スコアマッチングで行っているのは,\(\bm{x}\) に関するスコアマッチングではなく,ノイズ付加後の \(\tilde{\bm{x}}\) に関するスコアマッチングである.したがって,推定したスコアを用いた Langevin モンテカルロにもバイアスがある.Hyvärinen (2025) は,ノイズが正規分布の場合は, \[ \bm{x}_{t+1}=\tilde{\bm{x}}_{t}+\varepsilon\nabla_{\tilde{\bm{x}}}\log{p}(\tilde{\bm{x}}_{t})+\sqrt{2\varepsilon-\sigma^{2}}\bm{z}_{t},\quad\bm{z}_{t}\sim N(0,I) \]

とすることでバイアスを補正できることを示した.実際,両辺に \(\bm{u}_{t+1}\sim N(0,\sigma^{2}I)\) を足すと,\(\sqrt{2\varepsilon-\sigma^{2}}\bm{z}_{t}+\bm{u}_{t+1}\sim N(0,2\varepsilon I)\) なので, \[ \tilde{\bm{x}}_{t+1}=\tilde{\bm{x}}_{t}+\varepsilon\nabla_{\tilde{\bm{x}}}\log{p}(\tilde{\bm{x}}_{t})+\sqrt{2\varepsilon}\bm{z}_{t}',\quad\bm{z}_{t}'\sim N(0,I) \]

となり,ノイズ付加後の \(\tilde{\bm{x}}\) に関して正確な Langevin モンテカルロが得られる(なので,ノイズ成分 \(\bm{u}_{t+1}\) を引いた更新とも思える).特に \(\varepsilon=\sigma^{2}/2\) とした場合,更新式は \[ \bm{x}_{t+1}=\tilde{\bm{x}}_{t}+ \frac{\sigma^{2}}{2}\nabla_{\tilde{\bm{x}}}\log{p}(\tilde{\bm{x}}_{t}),\quad \tilde{\bm{x}}_{t}=\bm{x}_{t}+\bm{u}_{t},\quad\bm{u}_{t}\sim N(0,\sigma^{2}I) \]

で与えられる(ハーフ・デノイジング).

スコアマッチングの拡張

比マッチング

ここでは \(\mathfrak{X}=\{-1,1\}^{d}\) の場合を考える.離散集合上では勾配を考えることができないため,別のアイデアが必要である.比マッチングは Hyvärinen (2007) によって提案された,離散集合上へのスコアマッチングの拡張である.\(\bm{x}\in\mathfrak{X}\) に対し,\(\bm{x}_{-j}=(x_1,\dots,-x_j,\dots,x_d)\) とする.

Definition 2 \(\mathfrak{X}=\{-1,1\}^{d}\) 上の2つの確率分布 \(P,Q\) は,密度 \(p(\bm{x}),q(\bm{x})\) をもつとする. \(g(\infty)=0\) を満たす \(g:[0,\infty]\to[0,\infty)\) に対し, \[ \begin{align*} J_{RM}(P\,\|\,Q)&:=\frac{1}{2} E_P\left[ \sum_{j=1}^{d}\left\{ g\left(\frac{p(\bm{x})}{p(\bm{x}_{-j})}\right)-g\left(\frac{q(\bm{x})}{q(\bm{x}_{-j})}\right)\right\}^{2}\right.\\ &\hphantom{=}\left.+\sum_{j=1}^{d}\left\{ g\left(\frac{p(\bm{x}_{-j})}{p(\bm{x})}\right)-g\left(\frac{q(\bm{x}_{-j})}{q(\bm{x})}\right)\right\}^{2} \right] \end{align*} \]

と定義する.

\(g\) は0除算に対処するための関数で,Hyvärinen (2007) では,\(g(u)=1/(1+u)\) が提案されている.この形は,後でスコアを導出する際にも本質的な役割を果たす.

Proposition 3 \(g\) が単射で,\(p>0\) なら,\(J_{RM}(P\,\|\,Q)=0\iff P=Q\) が成立.

証明

Proof. \((\Leftarrow)\) は明らか.\((\Rightarrow)\) を示す.仮定より, \[ g\left(\frac{p(\bm{x})}{p(\bm{x}_{-j})}\right) =g\left(\frac{q(\bm{x})}{q(\bm{x}_{-j})}\right) \]

が任意の \(\bm{x}\in\mathfrak{X}\) に対して成立.\(g\) は単射なので, \[ \frac{p(\bm{x})}{p(\bm{x}_{-j})}=\frac{q(\bm{x})}{q(\bm{x}_{-j})} \]

が成立. \[ \frac{p(\bm{x})}{q(\bm{x})}=\frac{p(\bm{x}_{-j})}{q(\bm{x}_{-j})} \]

である.同様に,任意の \(k\) に対し, \[ \frac{p(\bm{x}_{-j})}{q(\bm{x}_{-j})}=\frac{p(\bm{x}_{-j,-k})}{q(\bm{x}_{-j,-k})} \]

が成立.これを繰り返すことで,ある \(c>0\) が存在し,任意の \(\bm{x}\in\mathfrak{X}\) に対し,\(p(\bm{x})/q(\bm{x})=c\) がいえる.\(p,q\) は密度なので,\(p=q\)

比マッチングのスコアを導出する.簡単のために,\(p(\bm{x})/p(\bm{x}_{-j})\)\(p/p_{-j}\) と略記する.また,\(g(u)=1/(1+u)\) とする.

途中式

\[ \begin{align*} 2J_{RM}(P\,\|\,Q) &=\sum\limits_{j=1}^{d}E_{P}\left[ \{g(p/p_{-j})-g(q/q_{-j})\}^{2}+\{g(p_{-j}/p)-g(q_{-j}/q)\}^{2} \right]\\ &=\sum\limits_{j=1}^{d}E_{P}\left[ g(p/p_{-j})^{2}+g(p_{-j}/p)^{2} \right]+ \sum\limits_{j=1}^{d}E_{P}\left[\\ g(q/q_{-j})^{2}+g(q_{-j}/q)^{2} \right]\\ &\hphantom{=}-2\sum\limits_{j=1}^{d}E_{P}\left[ g(p/p_{-j})g(q/q_{-j})+g(p_{-j}/p)g(q_{-j}/q) \right] \end{align*} \]

である.第3項を評価する. \[ \begin{align*} &E_{P}\left[ g(p/p_{-j})g(q/q_{-j})+g(p_{-j}/p)g(q_{-j}/q) \right]\\ &=E_{P}\left[ \frac{p_{-j}}{p+p_{-j}}g(q/q_{-j})\right]+ E_{P}\left[\frac{p}{p+p_{-j}}g(q_{-j}/q) \right]\\ &=\sum\limits_{x\in\mathfrak{X}}p_{-j} \frac{p}{p+p_{-j}}g(q_{-j}/q) +\sum\limits_{x\in\mathfrak{X}}p\frac{p}{p+p_{-j}}g(q_{-j}/q)\\ &=E_{P}[g(q_{-j}/q)] \end{align*} \]

また, \[ \begin{align*} g(q/q_{-j})^{2}+g(q_{-j}/q)^{2}-2g(q_{-j}/q) &=\frac{q_{-j}^{2}+q^{2}-2(q+q_{-j})q}{(q+q_{-j})^{2}}\\ &=\frac{-(q+q_{-j})^{2}+2q_{-j}^{2}}{(q+q_{-j})^{2}} \end{align*} \]

比マッチングのスコアは

\[ \begin{align*} 2J_{RM}(P\,\|\,Q) &=\sum\limits_{j=1}^{d}E_{P}\left[ g(p/p_{-j})^{2}+g(p_{-j}/p)^{2} \right]+ \sum\limits_{j=1}^{d}E_{P}\left[ -1+2\frac{q_{-j}^{2}}{(q+q_{-j})^{2}} \right]\\ &=2E_{P}\left[\sum\limits_{j=1}^{d}g(q/q_{-j})^{2}\right]+\text{Const} \end{align*} \]

である.

一般化 Hyvärinen スコア

境界項が無視できない場合でも適当な重み \(\bm{w}(\bm{x})\) を導入することでスコアマッチングを適用できる.

Definition 3 (一般化 Fisher ダイバージェンス) \(\mathfrak{X}\subset\mathbb{R}^{d}\) 上の2つの確率分布 \(P,Q\) は,微分可能な正の密度 \(p(\bm{x}),q(\bm{x})\) をもつとする. \(\bm{w}=(w_{1},\dots,w_{d})^{\top},w_{j}:\mathfrak{X}\to[0,\infty)\) とする. \[ J_\bm{w}(P\,\|\,Q):=\int_{\mathfrak{X}} p(\bm{x})\|\bm{w}(\bm{x})^{1/2}\odot\nabla_{\bm{x}}\log p(\bm{x})-\bm{w}(\bm{x})^{1/2}\odot\nabla_{\bm{x}}\log q(\bm{x})\|^2\d{\bm{x}} \]

を一般化 Fisher ダイバージェンスという.ここで,\(\odot\) は要素積.

Proposition 4 \(\mathfrak{X}\) が開連結なら,一般化 Fisher ダイバージェンスはダイバージェンスの性質を満たす.すなわち,

  • \(J_\bm{w}(P\,\|\,Q)\geq0\)
  • \(J_\bm{w}(P\,\|\,Q)=0\iff P=Q\)

Proof. Proposition 1 と同様.

スコアを導出する.

途中式

一般化 Fisher ダイバージェンスは以下のように展開できる. \[ \begin{align*} J_\bm{w}(P\,\|\,Q) &=\int_{\mathfrak{X}} p(\bm{x})\|\bm{w}(\bm{x})^{1/2}\odot\nabla_{\bm{x}}\log{p}(\bm{x})-\bm{w}(\bm{x})^{1/2}\odot\nabla_{\bm{x}}\log{q}(\bm{x})\|^{2}\d{\bm{x}}\\ &=\int_{\mathfrak{X}} p(\bm{x}) \left( \|\bm{w}(\bm{x})^{1/2}\odot\nabla_{\bm{x}}\log{p}(\bm{x})\|^{2} +\|\bm{w}(\bm{x})^{1/2}\odot\nabla_{\bm{x}}\log{q}(\bm{x})\|^{2}\right.\\ &\hphantom{=}\left.-2\langle\nabla_{\bm{x}}\log{p}(\bm{x}),\bm{w}(\bm{x})\odot\nabla_{\bm{x}}\log{q}(\bm{x})\rangle \right)\d{\bm{x}} \end{align*} \]

ここで,第3項は Green の定理(あるいは部分積分)が適用可能なら,

\[ \begin{align*} &\int_\mathfrak{X}p(\bm{x})\langle\nabla_{\bm{x}}\log{p}(\bm{x}),\bm{w}(\bm{x})\odot\nabla_{\bm{x}}\log{q}(\bm{x})\rangle\d{x}\\ &=\sum\limits_{j=1}^{d}\int_\mathfrak{X}p(\bm{x})w_{j}(\bm{x})\pd{\log{p}(\bm{x})}{x_{j}}\pd{\log{q}(\bm{x})}{x_{j}}\d{x}\\ &=\sum\limits_{j=1}^{d}\int_\mathfrak{X}w_{j}(\bm{x})\pd{p(\bm{x})}{x_{j}}\pd{\log{q}(\bm{x})}{x_{j}}\d{x}\\ &=-\sum\limits_{j=1}^{d}\int_\mathfrak{X}p(\bm{x})\pd{}{x_{j}}\left(w_{j}(\bm{x})\pd{\log{q}(\bm{x})}{x_{j}}\right)\d{x}+\sum\limits_{j=1}^{d}\int_{\partial{\mathfrak{X}}}p(\bm{x})w_{j}(\bm{x})\pd{\log{q}(\bm{x})}{x_{j}}\nu_{j}(\bm{x})\d{s} \end{align*} \]

と展開できる.

境界項が無視できるならば,

\[ \begin{align*} J_\bm{w}(P\,\|\,Q) &=E_{P}\left[\|\bm{w}(\bm{x})^{1/2}\odot\nabla_{\bm{x}}\log{q}\|^{2}+2\sum_{j=1}^{d}\partial_{j}(w_{j}\partial_{j}\log{q})\right]+E_{P}\Big[\|\bm{w}(\bm{x})^{1/2}\odot\nabla_{\bm{x}}\log{p}\|^{2}\Big] \end{align*} \]

となる.\(\|\bm{w}(\bm{x})^{1/2}\odot\nabla_{\bm{x}}\log{q}\|^{2}+2\sum_{j=1}^{d}\partial_{j}(w_{j}\partial_{j}\log{q})\) を一般化 Hyvärinen スコアという.

\(\bm{w}\)\(\partial{\mathfrak{X}}\)\(0\) になるように選ぶことで,境界項は無視できる. \(\mathfrak{X}=[0,\infty)^{d}\) において, Hyvärinen (2007) は,\(\bm{w}(\bm{x}):=\bm{x}\) を,Yu, Drton, and Shojaie (2019)\(w_{j}(\bm{x})=\min\{x_{j},c\}\) を提案している.Liu, Kanamori, and Williams (2022) は,有界な \(\mathfrak{X}\) において, \(w_{j}(\bm{x}):=\inf_{\bm{z}\in\partial\mathfrak{X}}\|\bm{x}-\bm{z}\|\) とすることを提案している5

数値例

開球で切断された多変量正規分布について,一般化スコアマッチング推定量とスコアマッチング推定量(最尤推定量)を比べてみる.ここで,重み \(w\)Liu, Kanamori, and Williams (2022) に従って, \(w_{j}(\bm{x})=1-\|\bm{x}\|\) とした.

Code
using Distributions
using LinearAlgebra
using Plots

using JuMP
using Ipopt

function ScoreMatching(μ_true,Σ_true,X)
    n = size(X,2)
    model = Model(Ipopt.Optimizer)
    set_optimizer_attribute(model, "print_level", 0) 

    @variable(model, Λ[i=1:2,j=1:2],start = inv(Σ_true)[i,j])
    @variable(model, μ[i=1:2],start= μ_true[i])

    @objective(model,Min,sum([
        (1-norm(X[:,i]))*(X[:,i]-μ)'*Λ*Λ*(X[:,i]-μ) + 2/norm(X[:,i]) * X[:,i]'*Λ*(X[:,i]-μ) - 2(1-norm(X[:,i]))*(Λ[1,1]+Λ[2,2])
    for i in 1:n]))

    @constraint(model, Λ[1,1]0)
    @constraint(model, Λ[2,2]0)
    @constraint(model, Λ[1,2]==Λ[2,1])
    @constraint(model, Λ[1,1]*Λ[2,2] - Λ[1,2]^20)

    optimize!(model)
    return value.(μ), inv(value.(Λ))
end

function Contour(x,μ,Σ)
    C = zeros(length(x),length(x))
    for i in keys(x)
        for j in keys(x)
            if norm([x[i],x[j]]) < 1
                z = [x[i],x[j]] - μ
                C[i,j] = exp(-1/2 * z'inv(Σ)*z)
            end
        end
    end
    return C'
end

function simu(μ_true,Σ_true)
    n = 1000
    N = 100n
    X = rand(MultivariateNormal(μ_true,Σ_true),N);
    ind = mapslices(x -> norm(x) < 1, X, dims=1) |> v -> vec(v)
    X = X[:,ind][:,1:n]

    μ_mle = sum(X,dims=2)/n |> vec
    Σ_mle = (X .- μ_mle) |> v -> v*v' /n
    μ_sm, Σ_sm = ScoreMatching(μ_true,Σ_true,X)
    
    x = -1.5:1e-2:1.5

    pl1 = plot(aspect_ratio=1.0,lim=[-1.5,1.5],legend=:bottomleft,tick=false,title="Truth")
    contourf!(x,x,Contour(x,μ_true,Σ_true),levels=30,cmap=:PuRd_3,colorbar=false)
    pl2 = plot(aspect_ratio=1.0,lim=[-1.5,1.5],legend=:bottomleft,tick=false,title="MLE")
    contourf!(x,x,Contour(x,μ_mle,Σ_mle),levels=30,cmap=:PuRd_3,colorbar=false)
    pl3 = plot(aspect_ratio=1.0,lim=[-1.5,1.5],legend=:bottomleft,tick=false,title="GSM")
    contourf!(x,x,Contour(x,μ_sm,Σ_sm),levels=30,cmap=:PuRd_3,colorbar=false)

    return plot(pl1,pl2,pl3,layout=(1,3),size=(750,300))
end
simu (generic function with 1 method)
Code
simu([-0.5,-0.8],[0.1 0;0 10])
Code
simu([0,0],[.2 0;0 .2])
Code
simu([1.2,1.2],[5 -1;-1 1])
Code
simu([0,0.5],[2 1;1 1])

球面上のスコアマッチング

間に合わず.Mardia, Kent, and Laha (2016) を紹介する予定だった. 球面には境界が存在しないため,境界項が自動的に0になる.

組成データのスコアマッチング

間に合わず.Scealy and Wood (2023) を紹介する予定だった. 組成データは,確率単体 \(\Delta^{d-1}=\{\bm{u}\in\mathbb{R}^{d}\mid\bm{u}^\top\bm{1}=1,\bm{u}\geq\bm{0}\}\) 上のデータである.\(\bm{u}^{1/2}\) を考えることで,\(\Delta^{d-1}\) ではなく,\(\mathbb{S}^{d-1}\cap[0,\infty]^{d}\) 上のスコアマッチングに帰着させるようだ.

無限次元指数型分布族のスコアマッチング

間に合わず.Sriperumbudur et al. (2017) を紹介する予定だった.

欠測がある場合のスコアマッチング

間に合わず.Uehara, Matsuda, and Kim (2020)Givens, Liu, and Reeve (2025) を紹介する予定だった.

その他

斉次ダイバージェンス

Takenouchi and Kanamori (2017) では, 正規化定数の評価が困難な離散分布のパラメータ推定法として,斉次ダイバージェンスの局所化に基づく方法を提案している.

\(\mathfrak{X}\) を離散集合として, \[ \begin{align*} \mathcal{M}&:=\{f:\mathfrak{X}\to[0,\infty)\mid \ang{f}<\infty,f\neq0\},\\ \mathcal{P}&:=\{f\in\mathcal{M}\mid \ang{f}=1\} \end{align*} \]

とする.ただし,\(\ang{f}:=\sum_{x\in\mathfrak{X}}f(x)\) である(計数測度に関する積分).

Definition 4 (斉次ダイバージェンス) 非負関数 \(D:\mathcal{M}\times\mathcal{M}\to[0,\infty)\) が斉次ダイバージェンスであるとは,

  • \(f\in\mathcal{M}\) に対し,\(D(f,f)=0\).
  • \(f,g\in\mathcal{M}\) と任意の \(c>0\) に対し,\(D(f,g)=D(f,cg)\)
  • \(D(f,g)=0\) ならば,ある \(c>0\) が存在し,\(f=cg\)

を満たすことである.

定義から,斉次ダイバージェンスに基づく推定は,正規化定数に依らない.

Proposition 5 (Hölder の不等式) \((S,\mathcal{A},\mu)\) は測度空間,\(f,g\)\(S\) 上の \(\mathbb{C}\) 値可測関数とする.任意の \(1/p+1/q=1\) を満たす \(p,q\in[1,\infty]\) に対し, \[ \|fg\|_{1}\leq\|f\|_{p}\|g\|_{q} \]

が成立.さらに,\(p,q\in(1,\infty),f\in L^{p}(\mu),g\in L^{q}(\mu)\) ならば,等号成立は \(|f|^{p}\)\(|g|^{q}\)\(L^{1}(\mu)\) 上で線形従属の場合に限る.すなわち,\(\alpha=\beta=0\) でないある \(\alpha,\beta\geq0\) が存在し,\(\alpha|f|^{p}=\beta|g|^{q}\)\(\mu\)-a.e. で成立する場合に限る.

\(f,g\in\mathcal{M},\gamma>0\) に対し,\(p=1+\gamma,q=(1+\gamma)/\gamma\) として,\(f,g^{\gamma}\) に Hölder の不等式を使うと, \[ \ang{f^{1+\gamma}}^{\frac{1}{1+\gamma}}\ang{g^{1+\gamma}}^{\frac{\gamma}{1+\gamma}}\geq\ang{fg^{\gamma}} \]

が得られる.等号成立は \(f\propto g\) に限る.この結果を用いて定義される斉次ダイバージェンスが擬球ダイバージェンス (psuedo-spherical divergence; PS divergence) である.

Definition 5 (擬球ダイバージェンス) \(\gamma>0\) に対し,擬球ダイバージェンス \(D_\gamma:\mathcal{M}\times\mathcal{M}\to[0,\infty)\)\[ D_\gamma(f,g):=\frac{1}{1+\gamma}\log{\ang{f^{1+\gamma}}}+\frac{\gamma}{1+\gamma}\log{\ang{g^{1+\gamma}}}-\log{\ang{fg^{\gamma}}} \]

で定義する.

擬球ダイバージェンスは \(\gamma\) ダイバージェンス としても知られている.

\(X_1,\dots,X_n\sim P\) とすると,経験分布 \(\mathbb{P}_{n}\) の密度 \(p_{n}(x)\)\[ p_{n}(x)=\sum_{x\in\mathcal{Z}_{n}}\frac{n_{x}}{n} \]

である.ただし,\(\mathcal{Z}_{n}:=\{x\in\mathfrak{X}\mid \exists{i}, X_{i}=x \},n_{x}=\sum_{i=1}^{n}1_{\{x\}}(X_i)\) である.

いま,密度 \[ q(x\mid\theta)=\frac{\tilde{q}_\theta(x)}{Z(\theta)},\quad Z(\theta)=\ang{\tilde{q}_\theta(x)} \]

をもつ分布 \(Q_{\theta}\) によって,\(P\) を推定することを考える.擬球ダイバージェンスによる経験推定を考えると,

\[ \begin{align*} D_\gamma(p_{n},\tilde{q}_{\theta}) &=\frac{\gamma}{1+\gamma}\log{\ang{\tilde{q}_\theta^{1+\gamma}}}-\log{\ang{p_{n}\tilde{q}_{\theta}^{\gamma}}}+\text{Const}\\ &=\frac{\gamma}{1+\gamma}\log{\sum_{x\in\mathfrak{X}}\tilde{q}_\theta(x)^{1+\gamma}}-\log{\sum_{x\in\mathcal{Z}_{n}}\frac{n_{x}}{n}\tilde{q}_{\theta}(x)^{\gamma}}+\text{Const} \end{align*} \]

の最小化を考えることになる.第1項の評価が課題となる.Takenouchi and Kanamori (2017) は局所化によって,この課題を解決した.

Definition 6 (局所擬球ダイバージェンス) \(\alpha,\alpha'\in\mathbb{R}\) は相異なる実数とする.\(p\in\mathcal{P}\) に対し,\(\mathcal{Z}=\{x\in\mathfrak{X}\mid p(x)>0\}\) とする.\(q\in\mathcal{M}\)\(\mathcal{Z}\) 上で正とする. \[ \begin{align*} LD_{\alpha,\alpha',\gamma}(p,q) &:=D_{\gamma}((p^{\alpha}q^{1-\alpha})^{\frac{1}{1+\gamma}}1_{\mathcal{Z}},(p^{\alpha'}q^{1-\alpha'})^{\frac{1}{1+\gamma}}1_{\mathcal{Z}})\\ &=\frac{1}{1+\gamma}\log\sum_{x\in\mathcal{Z}}p^{\alpha}q^{1-\alpha} +\frac{\gamma}{1+\gamma}\log\sum_{x\in\mathcal{Z}}p^{\alpha'}q^{1-\alpha'} -\log\sum_{x\in\mathcal{Z}}p^{\bar\alpha}q^{1-\bar\alpha} \end{align*} \]

を局所擬球ダイバージェンスという.ただし,\(\bar\alpha=(\alpha+\gamma\alpha')/(1+\gamma)\) である.

\(LD_{\alpha,\alpha',\gamma}(p,q)=0\) のとき,\(\mathcal{Z}\) 上で \(p^{\alpha}q^{1-\alpha}\propto p^{\alpha'}q^{1-\alpha'}\),すなわち \(\mathcal{Z}\) 上で \(p\propto q\) が成立する.

経験分布 \(p_n\) とモデル \(\tilde{q}_\theta\) の間の局所擬球ダイバージェンスを求めると, \[ \begin{align*} &LD_{\alpha,\alpha',\gamma}(p,\tilde{q}_{\theta})\\ &=\frac{1}{1+\gamma}\log\sum_{x\in\mathcal{Z}_{n}}\left(\frac{n_{x}}{n}\right)^{\alpha}\tilde{q}_{\theta}^{1-\alpha} +\frac{\gamma}{1+\gamma}\log\sum_{x\in\mathcal{Z}_{n}}\left(\frac{n_{x}}{n}\right)^{\alpha'}\tilde{q}_{\theta}^{1-\alpha'} -\log\sum_{x\in\mathcal{Z}_{n}}\left(\frac{n_{x}}{n}\right)^{\bar\alpha}\tilde{q}_{\theta}^{1-\bar\alpha} \end{align*} \]

どの項の和も最大 \(n\) 回で済むようになった.

Theorem 1 (Takenouchi and Kanamori (2017) Theorem 4) 任意の \(p\in\mathcal{P}\) と指数型分布族のモデル \(\tilde{q}_{\bm\theta}(x)=\exp(\bm{\theta}^\top\bm{t}(x))\) に対し,\(\bar\alpha=1\) ならば,\(LD_{\alpha,\alpha',\gamma}(p,\tilde{q}_{\bm\theta})\)\(\bm\theta\) に関して凸関数.

\(\bar\alpha=1\) となるように \(\gamma\) を選ぶと,

\[ LD_{\alpha,\alpha',\gamma}(p,q) =\frac{1-\alpha'}{\alpha-\alpha'}\log\sum_{x\in\mathcal{Z}}p^{\alpha}q^{1-\alpha} +\frac{\alpha-1}{\alpha-\alpha'}\log\sum_{x\in\mathcal{Z}}p^{\alpha'}q^{1-\alpha'} \]

となる.

Caution

局所擬球ダイバージェンスに基づく推定量は有効推定量になる.これは,\(\alpha\) ダイバージェンスと深く関係しているようだ,時間が足りずまとめられなかった.

\(\alpha,\alpha'\) の選び方

Takenouchi and Kanamori (2017) では,最尤推定量との関係に注目して \((\alpha,\alpha')\approx(1,0)\) の場合について,数値実験で詳しく調べている.モデル \(\tilde{q}_\theta(x)=\exp(\psi_\theta(x))\) を考える. 擬球ダイバージェンスに基づく推定量は \(\nabla_\theta LD_{\alpha,\alpha',\gamma}(p,\tilde{q}_\theta)=0\),すなわち

\[ \frac{\sum_{x\in\mathcal{Z}_{n}}(n_{x}/n)^{\alpha}\tilde{q}_\theta(x)^{1-\alpha}\nabla_\theta\psi_\theta(x)}{\sum_{x\in\mathcal{Z}_{n}}(n_{x}/n)^{\alpha}\tilde{q}_\theta(x)^{1-\alpha}} - \frac{\sum_{x\in\mathcal{Z}_{n}}(n_{x}/n)^{\alpha'}\tilde{q}_\theta(x)^{1-\alpha'}\nabla_\theta\psi_\theta(x)}{\sum_{x\in\mathcal{Z}_{n}}(n_{x}/n)^{\alpha'}\tilde{q}_\theta(x)^{1-\alpha'}}=0 \]

を満たす.\((\alpha,\alpha')=(1,0)\) を代入すると, \[ \sum_{x\in\mathcal{Z}_n}\frac{n_{x}}{n}\nabla_\theta\psi_\theta(x)- \frac{\sum_{x\in\mathcal{Z}_n}\tilde{q}_\theta(x)\nabla_\theta\psi_\theta(x)}{\sum_{x\in\mathcal{Z}_n}\tilde{q}_\theta(x)}=0 \]

となる.尤度方程式は \[ \sum_{x\in\mathcal{Z}_n}\frac{n_{x}}{n}\nabla_\theta\psi_\theta(x)- \frac{\sum_{x\in\mathcal{X}}\tilde{q}_\theta(x)\nabla_\theta\psi_\theta(x)}{\sum_{x\in\mathcal{X}}\tilde{q}_\theta(x)}=0 \] となる.

ノイズ対照推定

間に合わず.Gutmann and Hyvärinen (2010)Matsuda and Hyvärinen (2019) を紹介する予定だった.

スコアリングルールとしての特徴づけ

Parry, Dawid, and Lauritzen (2012)\(\mathfrak{X}=(a,b)\) における,正規化定数に依らないパラメータの推定の特徴づけを行った.

Definition 7 (スコアリングルール) 可測空間 \((\mathfrak{X},\mathscr{A})\) 上の確率測度の集合を \(\mathcal{P}\) とする.スコアリングルール \(S:\mathfrak{X}\times\mathcal{P}\to\mathbb{R}\) とは,任意の \(P\in\mathcal{P}\) に対し,\(S(\,\cdot\,,P):\mathfrak{X}\to\mathbb{R}\) が可測となる関数である.

Definition 8 (proper) スコアリングルール \(S\) が proper であるとは,任意の \(P\in\mathcal{P}\) に対し, \[ S(P,P)\leq S(P,Q)\quad \forall{Q}\in\mathcal{P} \] が成立することである. ただし,\(S:\mathcal{P}\times\mathcal{P}\to\mathbb{R}\) は, \[ S(P,Q):=\int_{\mathfrak{X}}S(x,Q)\d{P}(x) \] で定義される.特に,等号成立が \(P=Q\) に限る場合,\(S\) は strictly proper であるという.

Definition 9 (局所性) \(\mathfrak{X}\) 上の測度 \(\mu\) に対し,\(\mathcal{P}=\{P\mid p:=\d{P}/\d{\mu} \in C^{m}(\mathfrak{X}),p>0\}\) とする. スコアリングルール \(S\) が(非負整数 \(m\) に関して) \(m\) 次局所的 (\(m\)-local) であるとは,\(s\in C^{\infty}(\mathfrak{X}\times(0,\infty)\times\mathbb{R}^{m})\) が存在し, \[ S(x,Q)=s(x,q(x),q'(x),\dots,q^{(m)}(x)) \] と表せることである.ただし,\(q\)\(Q\) の密度.

Note

proper なスコアリングルールによって,エントロピー \(H(P):=S(P,P)\) とダイバージェンス \(D(P,Q):=S(P,Q)-H(P)\) が定義される.

Definition 10 (斉次性) \(m\) 次局所的なスコアリングルール \(S\) が(\(h\in\mathbb{Z}\) に関して) \(h\) 次斉次的 (\(h\)-homogeneous) であるとは,任意の \(x\in\mathfrak{X}\) に対し \(s(x,\,\cdot\,):(0,\infty)\times\mathbb{R}^{m}\to\mathbb{R}\)\(h\) 次斉次的であること,すなわち, \[ s(x,\lambda y_{0},\lambda y_{1},\dots,\lambda y_{m})=\lambda^{h}s(x,y_{0},y_{1},\dots,y_{m}) \] が任意の \(\lambda>0\)\((x,y_{0},y_{1},\dots,y_{m})\in\mathfrak{X}\times(0,\infty)\times\mathbb{R}^{m}\) に対して成立することである.

Note

斉次的なスコアリングルールに基づく統計的推測は,モデルの正規化定数に依らない.

Example 15 (対数スコア) \(S(x,Q):=\log{q}(x)\) とすると,\(H(P)=\int_{\mathfrak{X}}p(x)\log{p}(x)\d{x}\) は Shannon エントロピー.\(D(P,Q)=\int_{\mathfrak{X}}p(x)\log\frac{p(x)}{q(x)}\d{x}\) は Kullback-Leibler ダイバージェンス.対数スコアは1次局所的で,strictly proper なスコアリングルール.

Example 16 (Hyvärinen スコア) Hyvärinen スコアは \[ S(x,Q):=2\frac{q''(x)}{q(x)}-\left(\frac{q'(x)}{q(x)}\right)^{2} \] とかける.Hyvärinen スコアは2次局所的,0次斉次的なスコアリングルール.

proper なスコアリングルールの特徴づけのために,いくつかの微分作用素を定義する.

Definition 11 \(f(x,y_0,y_1,\dots,y_m)\in C^\infty(\mathfrak{X}\times(0,\infty)\times\mathbb{R}^m)\) に対し, \[ \begin{align*} Df&:=\pd{f}{x}+\sum_{j=0}^{m}y_{j+1}\pd{f}{y_{j}},\\ \Lambda f&:=\sum_{j=0}^{m}(-1)^{j}D^{j}\left[\pd{f}{y_{j}}\right],\\ Ef&:=\sum_{j=0}^{m}y_{j}\pd{f}{y_j} \end{align*} \]

と定義する.

\(h\) 次斉次的な \(\phi\in C^\infty(\mathfrak{X}\times(0,\infty)\times\mathbb{R}^d)\) 全体を \(\mathcal{F}_{m}^{h}\) として,\(\mathcal{F}^{h}:=\bigcup_{m=0}^{\infty}\mathcal{F}_{m}^{h}\) とする.\(\mathcal{F}^{h}\) は線形空間.このとき,\(\Lambda\)\(\mathcal{F}^{1}\) から \(\mathcal{F}^{0}\) への線形写像になっていることが,Lemma 1Lemma 2 から従う.

Lemma 1Lemma 2

Lemma 1 \(f\in\mathcal{F}_{m}^{h}\) とする. \[ \begin{align*} &\frac{\partial{f}}{\partial{x}}\in \mathcal{F}_{m}^{h} &\frac{\partial{f}}{\partial{y_{j}}}\in \mathcal{F}_{m}^{h-1},\quad j=0,\dots,m \end{align*} \]

Proof. 任意の \((x,y_{0},y_{1},\dots,y_{m})\in\mathfrak{X}\times(0,\infty)\times\mathbb{R}^{m}\) に対し,\(g(x,y_{0},y_{1},\dots,y_{m})=f(x,\lambda{y}_{0},\lambda{y}_{1},\dots,\lambda{y}_{m})\) とすると, \[ \begin{align*} &\hphantom{\iff}\frac{\partial{g}}{\partial{y}_{j}}(x,{y}_{0},{y}_{1},\dots,{y}_{m})=\frac{\partial{(\lambda^{h} f)}}{\partial{y}_{j}}(x,{y}_{0},{y}_{1},\dots,{y}_{m})\\ &\iff \lambda\frac{\partial{f}}{\partial{y}_{j}}(x,\lambda{y}_{0},\lambda{y}_{1},\dots,\lambda{y}_{m})=\lambda^{h}\frac{\partial{f}}{\partial{y}_{j}}(x,{y}_{0},{y}_{1},\dots,{y}_{m}) \end{align*} \] が成立.よって,\(\partial{f}/\partial{y}_{j}\)\(h-1\) 次斉次的.また, \[ \begin{align*} &\hphantom{\iff}\frac{\partial{g}}{\partial{x}}(x,{y}_{0},{y}_{1},\dots,{y}_{m})=\frac{\partial{(\lambda^{h} f)}}{\partial{x}}(x,{y}_{0},{y}_{1},\dots,{y}_{m})\\ &\iff \frac{\partial{f}}{\partial{x}}(x,\lambda{y}_{0},\lambda{y}_{1},\dots,\lambda{y}_{m})=\lambda^{h}\frac{\partial{f}}{\partial{x}}(x,{y}_{0},{y}_{1},\dots,{y}_{m}) \end{align*} \] より,\(\partial{f}/{\partial{x}}\)\(h\) 次斉次的.

Lemma 2 \(f\in \mathcal{F}_{m}^{h}\) とする. \[ Df\in \mathcal{F}_{m+1}^{h},\quad j=1,\dots,d. \]

Proof. 任意の \((x,y_{0},y_{1},\dots,y_{m+1})\in\mathfrak{X}\times(0,\infty)\times\mathbb{R}^{m+1}\) に対し,Lemma 1 より, \[ \begin{align*} Df(x,\lambda{y}_{0},\lambda{y}_{1},\dots,\lambda{y}_{m+1}) &=\frac{\partial{f}}{\partial{x}}(x,\lambda{y}_{0},\lambda{y}_{1},\dots,\lambda{y}_m)+\sum\limits_{j=1}^{m}\lambda{y}_{j+1}\frac{\partial{f}}{\partial{y}_{j}}(x,\lambda{y}_{0},\lambda{y}_{1},\dots,\lambda{y}_m)\\ &=\lambda^{h}\frac{\partial{f}}{\partial{x}}(x,y_{0},y_{1},\dots,y_{m})+\lambda^{h}\sum\limits_{j=1}^{m}{y}_{j+1}\frac{\partial{f}}{\partial{y}_{j}}(x,y_{0},y_{1},\dots,y_{m})\\ &=\lambda^{h}Df(x,y_{0},y_{1},\dots,y_{m+1}) \end{align*} \] が成立.

Theorem 2 (Euler の斉次関数定理) \[ f\in\mathcal{F}_{m}^{h}\iff Ef=hf. \]

証明

Proof. 任意の \((x,y_{0},y_{1},\dots,y_{m})\in\mathfrak{X}\times(0,\infty)\times\mathbb{R}^{m}\) に対し,\(g:(0,\infty)\to\mathbb{R}\)\(g(\lambda)=f(x,\lambda{y}_{0},\lambda{y}_{1},\dots,\lambda{y}_{m})-\lambda^{h}f(x,y_{0},y_{1},\dots,y_{m})\) で定める. \((\Rightarrow)\) 仮定より \(g=0\) で,特に \(g'=0\) である.

\((1)\enspace h=0\) の場合

\[ \begin{align*} g'(\lambda) &=\sum\limits_{j=0}^{m}y_{j}\frac{\partial{f}}{\partial{y_{j}}}(x,\lambda{y}_{0},\lambda{y}_{1},\dots,\lambda{y}_{m})=0 \end{align*} \] である.\(0=g'(1)=Ef(x,y_{0},y_{1},\dots,y_{m})\) が成立.

\((2)\enspace h\neq0\) の場合

\[ \begin{align*} g'(\lambda) &=\sum\limits_{j=0}^{m}y_{j}\frac{\partial{f}}{\partial{y_{j}}}(x,\lambda{y}_{0},\lambda{y}_{1},\dots,\lambda{y}_{m})-h\lambda^{h-1}f(x,y_{0},y_{1},\dots,y_{m})=0 \end{align*} \] である.\(0=g'(1)=Ef(x,y_{0},y_{1},\dots,y_{m})-hf(x,y_{0},y_{1},\dots,y_{m})\) が成立. \((\Leftarrow)\) 仮定より,\(\lambda{g}'(\lambda)=Ef(x,\lambda{y}_{0},\lambda{y}_{1},\dots,\lambda{y}_{m})-h\lambda^{h}{f}(x,y_{0},y_{1},\dots,y_{m})=hg(\lambda)\) が成立.

\((1)\enspace h=0\) の場合

\(\lambda g'(\lambda)=0\) より,\(g\) は定数.\(g(1)=0\) なので,\(g=0\) である.

\((2)\enspace h\neq0\) の場合

微分方程式 \(\lambda g'(\lambda)- hg(\lambda)=0\) を条件 \(g(1)=0\) の下で解く. \[ \begin{align*} &\hphantom{\iff}\frac{g'(\lambda)}{g(\lambda)}= \frac{h}{\lambda}\\ &\iff \int \frac{g'(\lambda)}{g(\lambda)}\,d\lambda=\int \frac{h}{\lambda}\,d\lambda+\text{Const}\\ &\iff \int \frac{g'(\lambda)}{g(\lambda)}\,d\lambda=\int \frac{h}{\lambda}\,d\lambda+\text{Const}\\ &\iff \int \frac{1}{g}\,dg=h\log{\lambda}+\text{Const}\\ &\iff \log{|g(\lambda)|}=\log{\lambda}^{h}+\text{Const}\\ &\iff g(\lambda)=\lambda^{h}\times\text{Const} \end{align*} \] \(g(1)=\text{Const}=0\) より,\(\text{Const}=0\) である.よって,\(g=0\) が成立.

\(D,\Lambda,E\) の意味

\(D\)\(\dd{S(x,Q)}{x}\) を考えるときに現れる.\(S\)\(m\) 次局所的とすると,微分の連鎖率より, \[ \begin{align*} \dd{S(x,Q)}{x} &=\pd{S(x,Q)}{x}+q'(x)\pd{S(x,Q)}{y_0}+\dots+q^{(m+1)}\pd{S(x,Q)}{y_m}\\ &=D[s](x,q(x),q'(x),\dots,q^{(m+1)}(x)). \end{align*} \]

\(\Lambda\) は汎関数 \(J(u)=\int_{a}^{b}F(x,u(x),\dots,u^{m}(x))\d{x}\) の最小化を変分法で解いたときに現れる Euler-Lagrange 方程式 \[ \pd{F}{u}-\dd{}{x}\pd{F}{u'}+\dots+(-1)^{m}\dd[m]{}{x}\dd{F}{u^{(m)}} =0 \]

\(\Lambda F=0\) で表現する.\(\Lambda\) が proper なスコアリングルールの特徴づけに関係する理由は,(直感的には)\(S(P,Q)\) の変分が \(P=Q\)\(0\) になることが,proper なスコアリングルールの必要条件になるから.

\(E\)\(\mathcal{F}^{h}\) から \(\mathcal{F}^{h}\) への線形写像になっていることが,Lemma 1 から従う.また,Theorem 2 より,\(\mathcal{F}_{m}^{h}\)\(E\) の固有値 \(h\) に関する固有空間である.

Theorem 3 (proper なスコアリングルールの生成) \((y_0,y_1,\dots,y_m)\mapsto\phi(x,y_0,y_1,\dots,y_m)\) が 各 \(x\in\mathfrak{X}\) で凹関数になるような,\(\phi\in\mathcal{F}^{1}\) に対し,\(s:=\Lambda\phi\) とすると,\(S\) はある \(\mathcal{P}\) に関して proper.

証明

Proof. \[ g_{j}(x):=\pd{\phi}{y_j}(x,q(x),q'(x),\dots,q^{(m)}(x)) \]

とする.このとき,

\[ \begin{align*} S(x,Q) &=\Lambda\phi(x,q(x),q'(x),\dots,q^{(m)}(x))\\ &=\sum\limits_{j=0}^{m}(-1)^{j}D^{j}\left[\frac{\partial\phi}{\partial{y}_{j}}\right](x,q(x),q'(x),\dots,q^{(m+j)}(x))\\ &=\sum\limits_{j=0}^{m}(-1)^{j}g_{j}^{(j)}(x) \end{align*} \]

と表せる.部分積分の公式から, \[ \begin{align*} S(P,Q) &=\int_{a}^{b}p(x)S(x,Q)\d{x}\\ &=\sum\limits_{j=0}^{m}(-1)^{j}\int_{a}^{b}p(x)g_{j}^{(j)}(x)\d{x}\\ &=\sum\limits_{j=0}^{m}\int_{a}^{b}p^{(j)}(x)g_{j}(x)\d{x} +\sum\limits_{j=1}^{m}(-1)^{j}\sum\limits_{l=0}^{j-1}(-1)^{l}\left[p^{(l)}(x)g_{j}^{(j-l-1)}(x)\right]_{a}^{b} \end{align*} \] が得られる.

(繰り返し)部分積分の公式

\(f,g\in C^{n}(\mathfrak{X})\) とする.このとき,\((fg^{(n-1)})'=f'g^{(n-1)}+fg^{(n)}\) である. また,\((f'g^{(n-2)})'=f''g^{(n-2)}+f'g^{(n-1)}\) より, \[ (fg^{(n-1)})'=(f'g^{(n-2)})'-f''g^{(n-2)}+fg^{(n)} \] が成立.これを繰り返して, \[ fg^{(n)}+(-1)^{n}f^{(n)}g=\sum\limits_{j=0}^{n-1}(-1)^{j}(f^{(j)}g^{(n-j-1)})' \] が得られる.よって, \[ \int_{a}^{b} f(x)g^{(n)}(x)\d{x}+(-1)^{n}\int_{a}^{b} f^{(n)}(x)g(x)\d{x} =\sum_{j=0}^{n}(-1)^{j}\left[f^{j}(x)g^{n-j-1}(x)\right]_{a}^{b} \]

\[ S_{0}^{\phi}(P,Q):=\sum\limits_{j=0}^{m}\int_{a}^{b}p^{(j)}(x)g_{j}(x)\d{x} \] として,\(S_{B}^{\phi}(P,Q):=S(P,Q)-S_{0}^{\phi}(P,Q)\) とする.Theorem 2 より, \[ \begin{align*} S_{0}^{\phi}(P,P) &=\sum\limits_{j=0}^{m}\int_{a}^{b}p^{(j)}(x)\frac{\partial\phi}{\partial{y}_{j}}(x,p(x),p'(x),\dots,p^{(m)}(x))\d{x}\\ &=\int_{a}^{b}(E\phi)(x,p(x),p'(x),\dots,p^{(m)}(x))\d{x}\\ &=\int_{a}^{b}\phi(x,p(x),p'(x),\dots,p^{(m)}(x))\d{x}\\ \end{align*} \] が成り立ち, \[ \begin{align*} D_{0}^{\phi}(P,Q) &:=S_{0}^{\phi}(P,Q)-S_{0}^{\phi}(P,P)\\ &=\int_{a}^{b}\sum\limits_{j=0}^{m}\left(p^{(j)}(x)-q^{(j)}(x)\right)\frac{\partial\phi}{\partial{y}_{j}}(x,q(x),q'(x),\dots,q^{(m)}(x))\d{x}\\ &\hphantom{=}-\int_{a}^{b}\left(\phi(x,p(x),p'(x),\dots,p^{(m)}(x))-\phi(x,q(x),q'(x),\dots,q^{(m)}(x))\right)\d{x} \end{align*} \] となる.\(\bm{y}'=(p(x),p'(x),\dots,p^{(m)}(x))^{\top},\bm{y}=(q(x),q'(x),\dots,q^{(m)}(x))^{\top}\) とおくと,\(D_{0}^{\phi}\) の被積分関数は \[ (\bm{y}'-\bm{y})^{\top}\nabla_{\bm{y}}\phi(x,\bm{y})-(\phi(x,\bm{y}')-\phi(x,\bm{y})) \] とかける.\(\phi\) は凹関数なので,これは非負.よって,境界項 \(S_{B}^{\phi}(P,Q)-S_{B}^{\phi}(P,P)\) が常に消えるような \(\mathcal{P}\) においては,\(D_{0}^{\phi}\) はダイバージェンス,すなわち \(S(x,Q)\) は proper である.

Example 17 (一般化 Hyvärinen スコアの生成) \(\phi(x,y_0,y_1)=-w(x)y_1^2/y_0\) は一般化 Hyvärinen スコアを生成する.

導出

\[ \begin{align*} \Lambda\phi &=\pd{\phi}{y_0}-D\left[\pd{\phi}{y_1}\right] \\ &=\pd{\phi}{y_0}-\frac{\partial^{2}\phi}{\partial{x}\partial{y_{1}}}-y_{1}\frac{\partial^{2}\phi}{\partial{y_{0}}\partial{y_{1}}}-y_{2}\pd[2]{\phi}{y_{1}}\\ &=w(x)\frac{y_{1}^{2}}{y_{0}^{2}}+w'(x)\frac{2y_{1}}{y_{0}}-y_{1}w(x)\frac{2y_{1}}{y_{0}^{2}}-y_{2}w(x)\frac{-2}{y_{0}}\\ &=w(x)\left(\frac{2y_{2}}{y_{0}}-\frac{y_{1}^{2}}{y_{0}^{2}}\right)+2w'(x)\frac{y_{1}}{y_{0}} \end{align*} \]

である.1次元の一般化 Hyvärinen スコアは

\[ \begin{align*} &\left(w(x)^{1/2}\dd{\log{q}(x)}{x}\right)^{2}+2\dd{}{x}\left(w(x)\dd{\log{q}(x)}{x}\right)\\ &=w(x)\frac{q'(x)^{2}}{q(x)^{2}}+2w(x)\dd[2]{\log{q}(x)}{x}+2w'(x)\frac{q'(x)}{q(x)}\\ &=w(x)\frac{q'(x)^{2}}{q(x)^{2}}+2w(x)\left(\frac{q''(x)}{q(x)}-\frac{q'(x)^{2}}{q(x)^{2}}\right)+2w'(x)\frac{q'(x)}{q(x)}\\ &=w(x)\left(\frac{q'(x)^{2}}{q(x)^{2}}-2\frac{q''(x)}{q(x)}\right)+2w'(x)\frac{q'(x)}{q(x)} \end{align*} \]

となり,一致することがわかる.

次に,一般化 Fisherダイバージェンス を導出する.

\[ \pd{\phi}{y_0}=w(x)\frac{y_1^2}{y_0^2}, \pd{\phi}{y_1}=-w(x)\frac{2y_1}{y_0} \]

である. \[ \begin{align*} S(P,Q) &=\sum\limits_{j=0}^{1}\int_{a}^{b}p^{(j)}(x)g_{j}(x)\,dx+\sum\limits_{j=1}^{1}(-1)^{j}\sum\limits_{l=0}^{j-1}(-1)^{l}\left[p^{(l)}(x)g_{j}^{(j-l-1)}(x)\right]_{a}^{b}\\ &=\int_{a}^{b}\left(p(x)g_{0}(x)+p'(x)g_{1}(x)\right)\,dx +\Big[-p(x)g_{1}(x)\Big]_{a}^{b}\\ &=\int_{a}^{b}w(x)\left(p(x)\frac{q'(x)^{2}}{q(x)^{2}}+p'(x)\frac{-2q'(x)}{q(x)}\right)\,dx +\left[-p(x)w(x)\frac{-2q'(x)}{q(x)}\right]_{a}^{b}\\ &=\int_{a}^{b}w(x)\left(p(x)\frac{q'(x)^{2}}{q(x)^{2}}-2p'(x)\frac{q'(x)}{q(x)}\right)\,dx +\left[2w(x)p(x)\frac{q'(x)}{q(x)}\right]_{a}^{b}\\\\ S(P,P) &=\int_{a}^{b}w(x)\left(-\frac{p'(x)^{2}}{p(x)}\right)\,dx+\Big[2w(x)p'(x)\Big]_{a}^{b} \\\\ S(P,Q)-S(P,P) &=\int_{a}^{b}p(x)w(x)\left(\frac{p'(x)}{p(x)}-\frac{q'(x)}{q(x)}\right)^{2}\,dx +\left[w(x)\left(2p(x)\frac{q'(x)}{q(x)}-2p'(x)\right)\right]_{a}^{b} \end{align*} \]

文献紹介

松田 (2025) は本資料の作成に際して特に参考した.

References

Givens, Josh, Song Liu, and Henry W. J. Reeve. 2025. “Score Matching With Missing Data.” arXiv. https://doi.org/10.48550/arXiv.2506.00557.
Gutmann, Michael, and Aapo Hyvärinen. 2010. “Noise-Contrastive Estimation: A New Estimation Principle for Unnormalized Statistical Models.” In Proceedings of the Thirteenth International Conference on Artificial Intelligence and Statistics, edited by Yee Whye Teh and Mike Titterington, 9:297–304. Proceedings of Machine Learning Research. Chia Laguna Resort, Sardinia, Italy: PMLR. https://proceedings.mlr.press/v9/gutmann10a.html.
Hyvärinen, Aapo. 2005. “Estimation of Non-Normalized Statistical Models by Score Matching.” Edited by Peter Dayan. Journal of Machine Learning Research 6: 695–709.
———. 2007. “Some Extensions of Score Matching.” Computational Statistics & Data Analysis 51 (5): 2499–2512. https://doi.org/10.1016/j.csda.2006.09.003.
———. 2025. “A Noise-Corrected Langevin Algorithm and Sampling by Half-Denoising.” arXiv. https://doi.org/10.48550/arXiv.2410.05837.
Klein, Natalie, Josue Orellana, Scott L. Brincat, Earl K. Miller, and Robert E. Kass. 2020. “Torus Graphs for Multivariate Phase Coupling Analysis.” The Annals of Applied Statistics 14 (2). https://doi.org/10.1214/19-AOAS1300.
Lin, Lina, Mathias Drton, and Ali Shojaie. 2016. “Estimation of High-Dimensional Graphical Models Using Regularized Score Matching.” Electronic Journal of Statistics 10 (1). https://doi.org/10.1214/16-EJS1126.
Liu, Song, Takafumi Kanamori, and Daniel J Williams. 2022. “Estimating Density Models with Truncation Boundaries Using Score Matching.”
Mardia, Kanti V., John T. Kent, and Arnab K. Laha. 2016. “Score Matching Estimators for Directional Distributions.” arXiv. https://doi.org/10.48550/arXiv.1604.08470.
Matsuda, Takeru, and Aapo Hyvärinen. 2019. “Estimation of Non-Normalized Mixture Models.” In Proceedings of the Twenty-Second International Conference on Artificial Intelligence and Statistics, edited by Kamalika Chaudhuri and Masashi Sugiyama, 89:2555–63. Proceedings of Machine Learning Research. PMLR. https://proceedings.mlr.press/v89/matsuda19a.html.
Parry, Matthew, A. Philip Dawid, and Steffen Lauritzen. 2012. “Proper Local Scoring Rules.” The Annals of Statistics 40 (1). https://doi.org/10.1214/12-AOS971.
Scealy, Janice L., and Andrew T. A. Wood. 2023. “Score Matching for Compositional Distributions.” Journal of the American Statistical Association 118 (543): 1811–23. https://doi.org/10.1080/01621459.2021.2016422.
Song, Yang, Sahaj Garg, Jiaxin Shi, and Stefano Ermon. 2019. “Sliced Score Matching: A Scalable Approach to Density and Score Estimation.” arXiv. https://doi.org/10.48550/arXiv.1905.07088.
Sriperumbudur, Bharath, Kenji Fukumizu, Arthur Gretton, Aapo Hyvärinen, and Revant Kumar. 2017. “Density Estimation in Infinite Dimensional Exponential Families.” arXiv. https://doi.org/10.48550/arXiv.1312.3516.
Takenouchi, Takashi, and Takafumi Kanamori. 2017. “Statistical Inference with Unnormalized Discrete Models and Localized Homogeneous Divergences.”
Uehara, Masatoshi, Takeru Matsuda, and Jae Kwang Kim. 2020. “Imputation Estimators for Unnormalized Models with Missing Data.” In Proceedings of the Twenty Third International Conference on Artificial Intelligence and Statistics, 831–41. PMLR. https://proceedings.mlr.press/v108/uehara20b.html.
Vincent, Pascal. 2011. “A Connection Between Score Matching and Denoising Autoencoders.” Neural Computation 23 (7): 1661–74. https://doi.org/10.1162/NECO_a_00142.
Yu, Shiqing, Mathias Drton, and Ali Shojaie. 2019. “Generalized Score Matching for Non-Negative Data.”
松田孟留. 2025. “非正規化モデルによる統計的推測—スコアマッチングとノイズ対照推定—.” 一般社団法人 日本統計学会. https://doi.org/10.11329/jjssj.54.177.

Footnotes

  1. このままだと,\(\mathfrak{X}=\mathbb{R}^{d}\) において,\(V=\{\bm{x}\mid\|\bm{x}\|=1\}\) としても切断分布と呼べる.このように次元が落ちる状況は切断分布とは言わない気がする.例えば,\(V\) は開集合としておけばこのような状況は含めずに済む.↩︎

  2. ここでは Hyvärinen (2005) よりやや一般化された状況で Hyvärinen スコアを導出する.↩︎

  3. ただし,\(\mathfrak{X},f,g\) はいくらかの条件を満たす必要がある.↩︎

  4. \(\mathfrak{X}=\mathbb{R}^{d}\) の場合,Green の定理をそのまま適用することはできないが,部分積分を行うことで同じ結果が得られる.↩︎

  5. この \(w_{j}\) は微分可能とは限らないが,弱微分は定義できる.↩︎